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The fee spin-1 Ising (BEG) model has a dense ferromagnetic (d/) ground 
state instead of the ferromagnetic ground state at low temperature region and 
exhibits the dense ferromagnetic (df) - ferromagnetic (F) phase transition for 
d = D/J = 2.9, k = K/J = -0.5, I = L/J = and /i = H/J = 0. The critical 
behavior of the dense ferromagnetic [df) - ferromagnetic (F) phase transition 
has been investigated using the cellular automaton cooling and heating algo- 
rithms. The universality class and the type of the dense ferromagnetic (df) - 
ferromagnetic {F) phase transition have been researched within the framework 
of the finite - size scaling, the power law relations and the probability distri- 
bution. The results show that the dense ferromagnetic- ferromagnetic phase 
transition is of the second order and the model shows universal second order 
Ising critical behavior at d = 2.9 parameter value through k = —0.5 line. 

PACS number(s): 05.10.-a, 05.50.-Hq, 64.60.-i 
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I . INTRODUCTION 

In recent years, some of the studies indicated that the spin-1 Ising model 
has a ground state ordered structure which is named the dense ferromagnetic 
( df ) [1 — 4]. In absence of the df ordered structure, the phase diagrams were 
formed for some phase transitions which has been assumed the weak first order 
transition instead of the second order transition. The presence of the df ordered 
structure can be caused to some changes on the global phase diagrams [4]. This 
case clarifies the differences among the results of the previous studies for the 
{kTc/J, d) phase diagram through the k line. [5 — 11]. For example, the phase 
diagram has exhibited a tricritical point {TCP) instead of a critical end point 
{CEP) for k = -0.5. While the MFA [7] and RG [8] studies exhibited a critical 
end point (CEP) for k=-0.5 on the BEG model global phase diagram, CA [5, 6], 
MCRG [9], TPCA [10] and CVM [11] studies show that there is the tricritical 
point (TCP). Through the k line, the df - F phase transition is very important 
for determining the type of the phase boundary. The purpose of this study is 
to define the df ordered phase and is to investigate the nature of df — F phase 
transition at d = 2.9 value through k = —0.5 line. This point creates the type 
of the special point of the ( kTc/J , d ) phase space for k = —0.5. Therefore, 
we have found the df - F phase transition worthy of investigation in depth. 
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Furthermore, the universaUty class of the df - F phase transition has not been 
investigated so far. The critical temperature and the statical critical exponents 
are estimated by analyzing the data within the framework of the finite - size 
scaling theory and the power law relations. 

The spin-1 Ising model, which is known as the generalized Blume-Emery- 
Griffiths (BEG) model, can be used to simulate many physical systems. The 
model firstly has been presented for describing phase separation and superfluid 
ordering in He mixtures [12]. The versions of the model have been applied to the 
physical systems such as the solid-liquid-gas systems [13], the multicomponent 
fluids [14], the microemulsions [15], the semiconductor alloys [16 — 18], He^-He* 
mixtures [12, 19] and the binary alloys [20]. 

The BEG model Hamiltonian is defined as 

<ij> <ij> <ij> i i 

(1) 

which is equivalent to the lattice gas Hamiltonian under some transformations 
[21 — 23]. {ij) denotes summation over all nearest-neighbor (nn) pairs of sites 
and Si = —1, 0, 1. The parameters J, K, L, D and h are bilinear, biquadratic, 
dipole-quadrupole interaction terms, the single-ion anisotropy constant and the 
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field term. The BEG model for fc > has been studied by mean field ap- 
proximation (MFA) [12 — 14], the transfer matrix method [24], series expan- 
sion method [25], the constant coupling approximation [26], the position-space 
renormalization method [27], cluster variation method (CVM) [1], linear-chain 
approximation [2] , Monte Carlo method (MC) [3] and Cellular Automaton (CA) 
[4 - 6, 28, 29]. 

In this paper, the fee BEG model ior d = D/J = 2.9, k = K/J = -0.5, 
I = Lj J = ^ and h = H/J = is simulated using cooling and heating algorithm 
improved from Creutz Cellular Automaton. In the previous papers, the Creutz 
cellular automaton (CCA) algorithm and its improved versions have been used 
successfully to study the properties of the critical behaviors of the Ising model 
Hamiltonians [4 — 6, 28 — 50]. The CCA algorithm, which was first introduced 
by Creutz [51], is a microcanonical algorithm interpolating between the conven- 
tional Monte Carlo and the molecular dynamics techniques. The Creutz cellular 
automaton (CCA) is faster than the conventional Monte Carlo method (MC). 
The CCA does not need high quality random numbers and it is a new and an 
alternative simulation method for physical systems. It has another advantage 
allowing the specific heat to be computed from internal energy fluctuations. 
Our previous studies showed that the heating and the cooling algorithms im- 
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proved from the Creutz Cellular Automaton algorithm are effective to study 
the phase space and the critical behavior of the Blume Emery Griffiths model 
[4-6,28-30,35,36]. 

II . RESULTS AND DISCUSSION 

The CA algorithm of spin-1 Ising model is a microcanonical algorithm. The 
total energy H, which is conserved, is given by 

H = Hi + Hk (2) 

where Hj is Ising energy which is given by equationl and Hk is kinetic energy. 

The kinetic energy Hk is an integer, equal to the change in the Ising spin 
energy for any spin flip and its value lie in the interval (0, m). m is equal to 
24J for d = D/J = 2.9, k = K/J = -0.5, £ = L/J = and h = H/J = on 
fee lattice. For a site to be updated, its spin is changed to one of the other two 
states with 1/2 probability. If this energy is transferable to or from the kinetic 
energy variable of the site, such that the total energy H is conserved, then this 
change is done and kinetic energy is appropriately changed. Otherwise the spin 
is not change [30, 34 — 36]. 

At the heating and the cooling algorithms, the simulation consist of two 
parts, the initialization procedure and the computation of the thermodynamic 
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quantities. The initial configuration for heating and coohng algorithms can be 
set in different shapes. In this study, the initial configurations are obtained at 
three different shapes for heating algorithm during 20.000 CA steps. Firstly, all 
the spins are up {S = +1) at the absolute zero temperature for both algorithms. 
The initial configuration of the heating algorithm has been obtained at low 
temperature ordered phase {df) adding kinetic energy which is equal to the 
maximum change in the Ising spin energy for the any spin flip to the spin 
system for set I and II. The another initial configuration (set III) has been 
obtained flipping 8% of the spins to 5 = state. The heating rate is realized by 
increasing of 8% in the kinetic energy {Hk) of 15% of the fee lattice for two sets 
at the computation of the thermodynamic quantities. At set III, the heating 
rate is realized by increasing of 8% in the kinetic energy (-fffc) of all site of the 
fee lattice. 

On the other hand, the initial configuration for the cooling algorithm is 
obtained adding energy to the 70% of the spin system for getting the disordered 
phase (P) at high temperature. During the cooling cycle, the cooling rate is 
realized by decreasing of 8% in the kinetic energy {Hf.) from 25% of the spin 
system. The initial configurations are run during the 20.000 Cellular Automaton 
time steps. Instead of resetting the starting configuration at each energy, it 
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was used the final configuration at a given energy as the starting point for 
the next at both heating and coohng algorithms. The computed values of the 
thermodynamic quantities (the order parameters (M, Q), the susceptibility (x), 
the Ising energy {Hj) and the specific heat (C)) are averages over the lattice 
and over the number of time steps (2.000.000 ) with discard of the first 100.000 
time steps during the cellular automaton develops [4 — 6, 28, 29]. 

They have been computed on the fee lattice with L = 8, 9, 10, 11 and 12 
(The total number of sites is N = 4L^) for periodic boundary conditions. The 
fee lattice was formed in a simple cubic (sc) lattice. (The total number of sites 
is N = AL^ = 6912 for L = 12 fee lattice, this total site number equals to 
L = 19 for the simple cubic lattice). The presented figures are set III (heating 
algorithm) results. 

The order parameters, the Ising energy, susceptibility and specific heat are 
calculated from 




(3) 



<ij> <ij> 



(4) 



X = N 



(m2) - (Mf 
kf 



(5) 
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where Uq is the ground state energy at kT/J = 0. 

The ferromagnetic (F) and the paramagnetic (P) phases can be determined 
with the average occupation of the states (-P±i,o)- As the projectors for the 
states S = +1, -1 and are P+i = ^S{S+1), P_i = ^S{S-1) and Pq = 1-5^, 
the average occupation of the states are (P+i) = ^{Q + M), (P-i) = ^{Q — M) 
and (Pq) = 1 — <3) respectively. With considering the average occupation of the 
states, another ferromagnetic phase can be determined as the dense ferromag- 
netic phase {df). 

Ferromagnetic (P): (P+i) ^ (P_i) ^ (Pq) ^ 0, (M ^ Q ^ 0) 
Dense ferromagnetic {df): (P_i) 0; (P+i) ^ (Pq) ^ 0, (M ^ Q 7^ 0) 
Paramagnetic (P): (P_i) = (P+i) ^ (Pq) ^ 0, (M = 0,Q ^ 0). 
II .1 Temperature Variations of Thermodynamic Quantities for the 
dF — F — P Phase Transitions 

The temperature variation of the order parameters (M, Q), the susceptibihty 
(x), the Ising energy (Hj) and the specific heat (C/) are illustrated in figure 1 for 
exhibiting the general aspect of the successive df — F — P phase transitions at 
d = 2.9 parameter value through k = —0.5 line. As it is seen in figure 1(a) and 



8 



figure 1(c), the order parameters and the Ising energy appear continuously for 
df —F and F — P phase transitions. Therefore both phase transitions are of the 
second order as functional behavior. The susceptibility (x) and the specific heat 
(C/, C)exhibit two peaks at Tci and Tc2 temperatures corresponding to df — F 
and F — P phase transitions (Figure 1(b), 1(d) and figure 2). It can be seen from 
figure 1(c), the functional change of the Ising energy from order to order {df 
—F) phase transition is different from the view of the order to disorder {F — P) 
phase transition. For df - F phase transition, the Ising energy difference (AC/) is 
greater than for the F - P phase transition. The estimated critical temperature 
from susceptibility and specific heat maxima is compatible with each other for 
F - P phase transition. But the critical temperature values for df - F phase 
transition are not compatible (Figure 1(b) and 2(a)). Therefore specific heat 
(C) has been recalculated for only spin-spin interaction energy (Figure 2(b)). 
The spin-spin interaction energy U is determined as 



The specific heat calculated from U can show a sharp peak for df — F phase 
transition. Because, the first sum (f/) in Ising energy [Ui) distinguish the S = 
+1 and —1 states. Indeed, the obtained infinite critical temperature (Tci(oo)= 




(7) 



<i3> 
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1.52 ±0.04) from the susceptibility (x) and specific heat (C) peak temperatures 
are compatible with each other. Tc2(oo) is obtained from the susceptibility (x) 
and specific heat (C/ and C) peak temperatures as 3.20 ± 0.02. 

The temperature variations of the (-P+i), {P-i) and (Pq) are given in figure 
3. The initial configuration for the heating algorithm is created as all spins are 
up {S = +1) at absolute zero temperature (T = 0). If the enough energy is 
added to the spin system, the 5 = begins to arise. The excitation energy of 
the single spin fiipping from 5 = +1 to 5 = is 3.1J while it is 24 J for the single 
spin flipping from S' = +ltoS'=— 1. Therefore, at low temperature region, 
the rising probability of the S = —1 state has to be lower than S = state. It 
can be seen in figure 3 that the spin system includes 5 = +1 and S = states 
predominantly. The value of (P+i) is about 1 and (Pq) is different from zero, 
while (P-i) appears almost zero indicating the df phase for T < Tci{L). At the 
same time,. So, M is almost equal to the Q (M = Q ^ 0). As {P-i) increases 
above Tci{L), the df ordered phase changes to the F ordered phase. At high 
temperature region, the ferromagnetic - the paramagnetic phase transition (F 
- P) occurs . (P+i) is equal to {P-i) and the system is in the paramagnetic 
(P) disordered phase above Tc2{L) temperature { M ^ Q ^ Q). Therefore the 
phase space is divided into three regions { df, F and P ). It is obvious that to 
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follow the temperature variation of {P-i) is a useful way to prove the existence 
of the df ordered phase. 

II .2 Probability Distribution of Order Parameter 

The another useful procedure to distinguish the phase transition type is to 
calculate the probability distributions of the order parameter (P(M)). In our 
study the probability distribution is calculated by 

P^i^) = (8) 
where Nm is the number of times that magnetization M appears, and Nccas is 
the total number of the cellular automaton steps. The histogram with 200 bins 
are used for plotting the probability distribution of the magnetization [24, 31]. 
The probability distribution of the order parameter (P(M)) near the phase 
transition temperature shows two peaks in the second phase transitions. 

The probability distributions of the order parameter (P(M)) are shown for 
different temperature values in figure 4. The peaks of the order parameter 
probability distribution exhibits minimum with increasing temperature at the 
low temperature region. This minimum corresponds to the second order df — F 
phase transition at the T^i{L = 12) = 1.479. Although the phase transition is 
of the second order, the probability distribution shows the single peak near the 
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phase transition temperature Tci{L). Because the system has S = +1 and 
spins below Tci{L) and the transition is from order (df) to order (F). However, 
the probability distributions in the F — P phase transition region exhibit the two 
peaks with the contribution of the S = —1 state near the T^^iL = 12) = 3.187. 
For T > Tc2, there is a single peak focused to M = indicating the disordered 
(P) phase. 

II. 3 Finite - Size Scaling Analyses and the Statical Critical Expo- 
nents 

The values of the statical critical exponents (u, j5, 7, a) are estimated within 
the framework of the finite - size scaling theory and the power laws. The infinite 
lattice critical temperature Tc (00) has been obtained from the susceptibility and 
the specific heat peak temperatures for the successive second order df —F — P 
phase transitions and from the intersection point of Binder cumulant curves {Ul) 
for the second order F —P phase transition [52]. The finite - size scaling relations 
of the Binder cumulant {Ul), the order parameter (M ), the susceptibility (x) 
and the specific heat (C) are given by 

Ul = G^ieL^'") (9) 
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(10) 



kTx = L'^^''Y°{sL^/'') (11) 

C = L"/''Z°{eL^/'') (12) 

For large x = eL^^", the finite lattice critical behaviors must be asymptoti- 
cally reproduced, that is, 

X°{x) oc Ax^ (13) 



Y°{x)(xBx-^ (14) 



Z°{x) oc Cx"" (15) 

According to the finite size scaling theory, the data for the finite - size 
lattices of the thermodynamic quantities should lie on a single curve for the 
temperatures both below and above Tc(oo) with universal critical exponents. 
The critical exponents 13, 7 and a have been obtained from the Log-Log plots 
of the asymptomatic functions. 
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The temperature variation of the Binder cumulant is shown for the different 
lattice sizes in figure 5. The Binder cumulant curves intersect at the Tq^^oo) = 
3.20±0.02 corresponding to the F — P phase transition (Figure 5(a)). This value 
is compatible with r^2(oo) which is extrapolated according to the finite size 
scaling theory from the susceptibility and the specific heat peak temperatures 
{Tc{L)) of the finite lattices, respectively. 

Tc(i) = Tc(oo) + oL-^/" (16) 

It can be seen in the inset of the figure 5(a) that there is no intersection at the 
Binder cumulant for the data of the df — F phase transition region. However, 
the Binder cumulant curves exhibit a plateau near the infinite lattice critical 
temperature {Tci{ooj) which is obtained from the susceptibility (x) and the 
specific heat (C) peak temperatures as T^j(oo) = 1.52 ± 0.04. In figure 5(b), 
the scaling data of the Binder cumulants are shown for the second order phase 
transition from F to P. Near the Tq2{oo), Binder cumulant curves have been 
scaled well for Tc=Tc2{oo) with u = 0.64. It can be seen in the inset of the 
figure 5(b) that the data corresponding to df ordered phase could not be scaled 
with the (°°) critical temperature value. However, the scaling data of the 
binder cumulant corresponding to the df ordered phase lie on a single curve for 
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e =(T - T^i(oo))/T^i(oo) at T < Tci{oo) region using u = 0.64 (Figure 5(c)) 
and the finite size scaling relations validate for the df ordered region. 

In figure 6(a), the scaling data of the order parameter is illustrated at the suc- 
cessive df — F — P phase transitions for L = 8, 9, 10, 11 and 12 at = T^2(°o)- 
The order parameter data lie on the two diflferent curves with slope=/3/z/ = 0.31 
and (3' /v = —0.55 for the temperatures both below and above Tq2{oo) respec- 
tively except for df ordered region with (3 = 0.31 and u = 0.64. As it is seen 
in the inset of the figure 6(a), the data of df — F phase transition region have 
not been scaled with Tc = Tq2{oo). However, the data corresponding to the 
df — F phase transition region have been scaled well with Tc = Tqi{oo) for 
T < Tci(oo) using /3 = 0.31 and ly = 0.64 in figure 6(b). 

The scaling data of the susceptibility have been shown in figure 7 with the 
straight lines describing the theoretically predicted behavior for large x (Equa- 
tion 8). The susceptibility data for the temperatures both below and above 
Tf-,2{oo) agrees with the asymptotic form except for the df - F phase transition 
region and so with the Tc = T^2(°°)' 7 = 7' = 1-25 and = 0.64 in figure 7(a) 
and 7(c). However, the data of df ordered phase (T < Tci) have been scaled 
with the Tc = T^i(oo) using 7 = 1.25 and v = 0.64 in figure 7(b). 

The finite size scaling data of the singular portion of the specific heat (C) 
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have been exhibited in figure 8. The data of F — P phase transition region 
of C/ have been scaled well both below and above TQ^ioo) using a = 0.12, 
f = 0.64 and the correction terms, b~ = —70 and 6+ = —8 (Figure 8(a) and 
8(b)). Although the data of the df phase region could not been scaled with Tc = 
Tq2{oo) in figure 8(a), the data lie on single curves with slope=— a/i/ = —0.12 
at the both side of the T^j(oo), using a = 0.12, v = 0.64 and the correction 
terms, b~ = —1 and 6+ = —0.3 in figure 8(c) and (d). On the other hand, 
the specific heat (Cj) data calculated from Uj scales well at df - F- P phase 
transitions using a = 0.12, u = 0.64 for T < r^^((X))and T < Tg^{oo) and 
T>T^2^(oo). 

The M, X and C data have been analyzed within the framework of the finite 
size scaling theory for the successive df — F — P phase transitions. The estimated 
values of the statical critical exponents are in good agreement with the universal 
values {a = 0.12, l3 = 0.31, 7 = 1.25, ly = 0.64) for the df - F and the F - P 
phase transitions. 

II. 4 Power Law Relations and the Infinite Lattice Statical Critical 
Exponents 

On the other hand, the critical exponent values for df — F phase transition 
can be obtained using the following power law relations [53]. 
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M{L) = sl^^^^ 



(17) 



X{L) = 



(18) 



C{L) = e 



-a(L) 



(19) 



where e ={T - Tc{L))/Tc{L). The finite lattice critical exponents /3(L), /3'(L), 
7(i), 7'(i), Oi{L) and Q!'(L) of the order parameter (M), susceptibility (x) and 
the specific heat (C) quantities are obtained from the slope of the log-log plot 
of the power laws relations for each finite lattices in the interval 0.05 < e < 0.2. 
The infinite lattice critical exponents are obtained using linear extrapolation and 
their values are given in Table I. The estimated values for cooling algorithm and 
three simulation sets of heating algorithm are in good agreement with the finite 
size scaling critical exponent estimations and the universal values for 3d Ising 
model {/3 = 0.31, 7 = 1.25, a = a' = 0.12 and u = 0.64). 

Table 1. The estimated values of the infinite lattice critical exponents and the 
critical temperatures (a, /?, 7 and Tq'^ (po)) using lineax extrapolation. 
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df - F phase transition 




Heating 


Cooling 




Set I 


Set II 


Set III 


Average of sets 






1.50 ±0.02 


1.50 ±0.03 


1.52 ± 0.03 


1.51 ±0.03 


1.52 ±0.04 








1.52 ±0.01 




1.52 ±0.01 


/3(T < Tci) 


0.31 ±0.01 


0.31 ±0.01 


0.30 ±0.01 


0.31 ±0.01 


0.31 ±0.01 


7(T < Tci) 


1.25 ±0.01 


1.23 ±0.03 


1.23 ±0.02 


1.24 ±0.03 


1.23 ±0.02 


a(T < Tci) 


0.12 ±0.01 


0.12 ±0.01 


0.11 ±0.01 


0.12 ±0.01 


0.12 ±0.01 


a'(T > Tci) 


0.12 ±0.01 


0.12 ±0.01 


0.12 ±0.01 


0.12 ±0.01 


0.12 ±0.01 



III. SUMMARY 

The {kTc / J, d) phase diagrams is illustrated in figure 9 for the presence of 
df order. The type of special point is determined by the d = 2.9 parameter. The 
calculations show that model exhibits the phase transition from order to order 
for d = 2.9 and the first order phase transition from order to disorder in the 
3 < d < 4 parameter region. If the model doesn't exhibit the df ordered phase 
instead of the F ordered phase in the low temperature region, the continuous 
phase transition from F order to F order {F — F) is considered as the weak 
first order. This constitutes the part of the first order phase transition line 
which creates the critical end point {CEP) [7,8]. However, CA results show 
that, the model has a df ordered phase for the parameters in the 2.9 < < 4.0 
region. The spin system contains S = +1 and S = states. As a result of this, 
the order parameters M and Q are almost equal each other (M = Q ^ 0) at 
low temperatures. With increasing temperature, the dense ferromagnetic (df) 
ordered phase changes continuously to ferromagnetic {F) ordered phase at the 
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d = 2.9, k = —0.5 with the enough contribution of S = —1 state. Therefore, near 
the d = 2.9, the first order phase transition fine have been changed to second 
order phase transition line, and there occurs the tricritical point TCP (Figure 
9). In order to determine the universahty class of the successive df — F — P 
second order phase transitions, the static critical exponents ( a, /3, 7 and v) are 
estimated within the framework of the finite - size scaling theory. The estimated 
values of the critical exponents ( a = 0.12, /3 = 0.31, 7 = 1.25 and v = 0.64) 
near the Tci and Tc2 temperatures are in good agreement with the theoretical 
values for three sets. The df — F phase transition is analyzed with the power 
laws for comparing with the critical exponent values estimated from the finite 
- size scaling theory. The obtained values are in compatible with the finite - 
size scaling analyze results and the universal values for the 3d Ising model. The 
obtained results have shown that the df — F phase transition is of the second 
order and it is compatible with the universal Ising critical behavior for d = 2.9 
parameter value through k = —0.5 line. As a result of this, the definition of the 
df phase changes the phase transition type and the special point type in the 
phase space for the BEG model. This result will lead to reexamine the structure 
of phase spaces. 
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Figure Captions 

Figure 1. For {d = 2.9, k = —0.5), the temperature dependence of (a) the 
order parameters (M, Q), (b) the susceptibility (x), (c) the Ising energy {Hi) 
and (d) the specific heat {Ci/k). 

Figure 2. For [d = 2.9, k = —0.5), the specific heat (C/k) calculated firom 

U. 

Figure 3. For {d = 2.9, k = —0.5), the temperature dependences of the 
(P). {P+i), {P-i) and (Pq) correspond to the 5 = +1, —1 and spin states, 
respectively. 

Figure 4. The probability distribution of the M for {d = —0.5, k = 0.9) on 
L = 12. 

Figure 5. For {d = 2.9, k = —0.5), (a) the temperature dependence of the 
Binder cumulant {Ul), (b) the finite - size scaling of the Binder cumulant near 

23 



the df — F — P phase transition with Tq2{oo), (c) the finite - size scaling of 
the Binder cumulant near the df — F phase transition with Tq-^^{oo). 

Figure 6. For {d = 2.9, k = —0.5), the finite - size scaUng plots of (a) the 
order parameter with T^2(oo)! (b) the order parameter with T^-^(oo) near the 
df — F phase transition for T < T^-^^{oo). 

Figure 7. For {d = 2.9, k = —0.5), the finite - size scaling plots of the 
susceptibility (a) with e —{T — TQ^ioo)) /Tq2{oo) for T < T^j^i^)^ (b) near the 
df - F phase transition with s ={T - T^i(oo))/T^i(oo) for T < T^iioo), (c) 
with £ =(T - Tc2{oo))/T for T > Tc2{oo). 

Figure 8. For [d = 2.9, k = —0.5), the finite - size scaling plots of the specific 
heat (a) for T < T^aloo) with Tc2(oo), (b) for T > Tc2{oo) with T^aCoo), (c) 
for T < T^i(oo) with T^i(oo), (d) for T > T^i(oo) with T^i(oo). 

Figure 9. The phase diagram for k = —0.5. The phase space with df phase 
contains a TCP at = 2.9. 
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